Option Strict Off
Imports System.Xml
Imports System.Windows.Forms
Imports MatrixDataStructures
Imports twoDMatrixLibrary


Public Class frmDTIParamertersSet
    Dim NumberOfGradients As Integer = 6
    Dim GradientStrength, SmallDelta, BigDelta As Single
    Dim LoadedDTIDirections As twoDMatrixLibrary.Matrix
    Public Nx, Ny, Nz As Integer
    Dim SubtractedMatrix As Matrix4Ddouble
    Dim DiffusionMatrix As Matrix4Ddouble
    Dim EigenVector1, EigenVector2, EigenVector3 As Matrix4DdoubleS
    Dim EigenValue1, EigenValue2, EigenValue3 As Matrix3DSingle
    Dim FAMatrix As Matrix3DSingle
    Dim RawDTIData As Matrix4DInt16

    Public Event FACalculationComplete As EventHandler

#Region "Old Computation"
    Public Sub StartDTIAnalysis()
        Dim i As Integer
        Dim logmatrix As New Matrix4Ddouble(Me.RawDTIData.i, Me.RawDTIData.j, Me.RawDTIData.k, Me.RawDTIData.l)
        'Replace all zeros by 1
        logmatrix.Replace(0, 1)
        logmatrix.data = Me.RawDTIData.Log(False)
        Me.SubtractedMatrix = New Matrix4Ddouble(Me.RawDTIData.i, Me.RawDTIData.j, Me.RawDTIData.k, Me.RawDTIData.l)
        For i = 1 To Me.NumberOfGradients
            Me.SubtractedMatrix = load_slices_into_volume(subtract_volumes(load_slices_from_volume(logmatrix, 0), load_slices_from_volume(logmatrix, i)), Me.SubtractedMatrix, i - 1)
        Next

    End Sub
    Public Function load_slices_into_volume(ByVal slices_matrix As Matrix3DSingle, ByVal volumes_matrix As Matrix4Ddouble, ByVal volume_no As Integer) As Matrix4Ddouble
        Dim a, b, c As Integer
        For a = 0 To slices_matrix.x - 1
            For b = 0 To slices_matrix.y - 1
                For c = 0 To slices_matrix.z - 1
                    volumes_matrix.data(a, b, c, volume_no) = slices_matrix.data(a, b, c)
                Next
            Next
        Next
        Return volumes_matrix
    End Function
    Public Function load_slices_from_volume(ByVal data As Matrix4Ddouble, ByVal volume_no As Integer) As Matrix3DSingle
        Dim ans As New Matrix3DSingle(data.i, data.j, data.k)
        Dim i, j, k As Integer
        For i = 0 To data.i - 1
            For j = 0 To data.j - 1
                For k = 0 To data.k - 1
                    ans.data(i, j, k) = CSng(data.data(i, j, k, volume_no))
                Next
            Next
        Next
        Return ans
    End Function
    'Public Function log(ByVal twoDMatrixLibrary.Matrix As twoDMatrixLibrary.Matrix4DInt16) As twoDMatrixLibrary.Matrix4Ddouble
    '    Dim ans As New twoDMatrixLibrary.Matrix4Ddouble(matrix.i, twoDMatrixLibrary.Matrix.j, twoDMatrixLibrary.Matrix.k, twoDMatrixLibrary.Matrix.l)
    '    Dim i, j, k, l As Integer
    '    For i = 0 To twoDMatrixLibrary.Matrix.i - 1
    '        For j = 0 To twoDMatrixLibrary.Matrix.j - 1
    '            For k = 0 To twoDMatrixLibrary.Matrix.k - 1
    '                For l = 0 To twoDMatrixLibrary.Matrix.l - 1
    '                    ans.data(i, j, k, l) = Math.Log(matrix.data(i, j, k, l))
    '                Next
    '            Next
    '        Next
    '    Next
    '    Return ans
    'End Function
    Public Function subtract_volumes(ByVal Matrix1 As MatrixDataStructures.Matrix3DSingle, ByVal Matrixn As MatrixDataStructures.Matrix3DSingle) As MatrixDataStructures.Matrix3DSingle
        Dim ans As New MatrixDataStructures.Matrix3DSingle(Matrix1.x, Matrix1.y, Matrix1.z)
        'do twoDMatrixLibrary.Matrix 1 - twoDMatrixLibrary.Matrix n
        Dim i, j, k As Integer
        For i = 0 To Matrix1.x - 1
            For j = 0 To Matrix1.y - 1
                For k = 0 To Matrix1.z - 1
                    ans.data(i, j, k) = Matrix1.data(i, j, k) - Matrixn.data(i, j, k)
                Next
            Next
        Next
        Return ans
    End Function

    'Public Sub dmatrix_compute()
    '    Dim x As New twoDMatrixLibrary.Matrix2Dsingle(Me.NumberOfGradients, 6)
    '    Dim i As Integer
    '    For i = 0 To Me.NumberOfGradients - 1
    '        x.data(i, 0) = Me.LoadedDTIDirections(i, 0) * Me.LoadedDTIDirections(i, 0)
    '        x.data(i, 1) = Me.LoadedDTIDirections(i, 1) * Me.LoadedDTIDirections(i, 1)
    '        x.data(i, 2) = Me.LoadedDTIDirections(i, 2) * Me.LoadedDTIDirections(i, 2)
    '        x.data(i, 3) = 2 * Me.LoadedDTIDirections(i, 0) * Me.LoadedDTIDirections(i, 1)
    '        x.data(i, 4) = 2 * Me.LoadedDTIDirections(i, 1) * Me.LoadedDTIDirections(i, 2)
    '        x.data(i, 5) = 2 * Me.LoadedDTIDirections(i, 0) * Me.LoadedDTIDirections(i, 2)
    '    Next

    '    'a=x
    '    Dim a_t As twoDMatrixLibrary.Matrix2Dsingle = transpose(x)
    '    Dim ls As twoDMatrixLibrary.Matrix2Dsingle = transpose(x)
    '    ls = multiply(ls, x)
    '    gaussjinvert(ls)

    '    Dim l As twoDMatrixLibrary.Matrix2Dsingle = multiply(ls, a_t)
    '    ls.dispose()
    '    a_t.dispose()
    '    MsgBox(l.x)
    '    MsgBox(l.y)
    '    'reuse raw data matrices in DTI object as new diffusion twoDMatrixLibrary.Matrix
    '    Me.DiffusionMatrix = New twoDMatrixLibrary.Matrix4Ddouble(Me.Nx, Me.Ny, Me.Nz, 6)

    '    Dim res As New twoDMatrixLibrary.Matrix3DSingle(Nx, Ny, Nz)
    '    Dim trd As New twoDMatrixLibrary.Matrix3DSingle(Nx, Ny, Nz)

    '    Dim j, k As Integer
    '    For i = 0 To Me.Nx - 1
    '        For j = 0 To Me.Ny - 1
    '            For k = 0 To Me.Nz - 1
    '                Dim temp As New twoDMatrixLibrary.Matrix2Dsingle(Me.NumberOfGradients, 1)

    '                Dim n As Integer
    '                For n = 0 To Me.NumberOfGradients - 1
    '                    temp.data(n, 0) = Me.SubtractedMatrix.data(i, j, k, n)
    '                Next

    '                Dim dtemp As twoDMatrixLibrary.Matrix2Dsingle = multiply(l, temp)
    '                ls.dispose()
    '                Dim dif As twoDMatrixLibrary.Matrix2Dsingle = multiply(x, dtemp)
    '                dif = subtract(temp, dif)

    '                res.data(i, j, k) = ((multiply(transpose(dif), dif).data(0, 0) ^ 0.5)) / ((Me.NumberOfGradients - 5) ^ 0.5)
    '                trd.data(i, j, k) = dtemp.data(0, 0) + dtemp.data(1, 0) + dtemp.data(2, 0)

    '                If trd.data(i, j, k) / res.data(i, j, k) > 0 Then 'check for threshold
    '                    Me.DiffusionMatrix.data(i, j, k, 0) = dtemp.data(0, 0)
    '                    Me.DiffusionMatrix.data(i, j, k, 1) = dtemp.data(1, 0)
    '                    Me.DiffusionMatrix.data(i, j, k, 2) = dtemp.data(2, 0)
    '                    Me.DiffusionMatrix.data(i, j, k, 3) = dtemp.data(3, 0)
    '                    Me.DiffusionMatrix.data(i, j, k, 4) = dtemp.data(4, 0)
    '                    Me.DiffusionMatrix.data(i, j, k, 5) = dtemp.data(5, 0)
    '                Else

    '                End If

    '            Next
    '        Next
    '    Next

    'End Sub
    'Function transpose(ByVal data As System.Drawing.Drawing2D.Matrix) As System.Drawing.Drawing2D.Matrix
    '    Dim ans As New System.Drawing.Drawing2D.Matrix
    '    Dim i, j As Integer
    '    Dim indice(1), tindice(1) As Integer
    '    For i = 0 To data.Elements.GetLongLength(0) - 1
    '        For j = 0 To data.Elements.GetLongLength(1) - 1
    '            indice(0) = i
    '            indice(1) = j
    '            tindice(0) = j
    '            tindice(1) = i
    '            ans.Elements.SetValue(data.Elements.GetValue(indice), tindice)
    '        Next
    '    Next
    '    Return ans
    'End Function
    'Function transpose(ByVal data As twoDMatrixLibrary.Matrix2Dsingle) As twoDMatrixLibrary.Matrix2Dsingle
    '    Dim ans As New twoDMatrixLibrary.Matrix2Dsingle(data.y, data.x)
    '    Dim i, j As Integer
    '    For i = 0 To data.x - 1
    '        For j = 0 To data.y - 1
    '            ans.data(j, i) = data.data(i, j)
    '        Next
    '    Next
    '    Return ans
    'End Function
    'Sub gaussjinvert(ByVal a As twoDMatrixLibrary.Matrix2Dsingle)
    '    'copied from somewhere
    '    Dim mayor, tmp, pivinv As Double
    '    Dim n As Integer = a.x
    '    Dim indco(n - 1), indfi(n - 1), piv(n - 1) As Integer
    '    Dim i, j, k, t, h, col, fil As Integer

    '    For j = 0 To n - 1
    '        piv(j) = 0
    '    Next j
    '    For i = 0 To n - 1
    '        mayor = 0.0
    '        For j = 0 To n - 1
    '            If piv(j) <> 1 Then
    '                For k = 0 To n - 1
    '                    If piv(k) = 0 Then
    '                        If Math.Abs(a.data(j, k)) >= mayor Then
    '                            mayor = Math.Abs(a.data(j, k))
    '                            fil = j
    '                            col = k
    '                        End If
    '                    ElseIf piv(k) > 1 Then
    '                        MsgBox("Singular twoDMatrixLibrary.Matrix")
    '                        Exit Sub
    '                    End If
    '                Next
    '            End If
    '        Next
    '        piv(col) = piv(col) + 1
    '        If fil <> col Then
    '            For t = 0 To n - 1
    '                tmp = a.data(fil, t)
    '                a.data(fil, t) = a.data(col, t)
    '                a.data(col, t) = tmp
    '            Next
    '            '            tmp = x(fil)
    '            '            x(fil) = x(col)
    '            '            x(col) = tmp
    '        End If
    '        indfi(i) = fil
    '        indco(i) = col
    '        If a.data(col, col) = 0 Then
    '            MsgBox("Singular twoDMatrixLibrary.Matrix")
    '            Exit Sub
    '        End If
    '        pivinv = 1.0 / a.data(col, col)
    '        a.data(col, col) = 1.0
    '        For t = 0 To n - 1
    '            a.data(col, t) = a.data(col, t) * pivinv
    '        Next
    '        'x(col) = x(col) * pivinv
    '        For h = 0 To n - 1
    '            If h <> col Then
    '                tmp = a.data(h, col)
    '                a.data(h, col) = 0.0
    '                For t = 0 To n - 1
    '                    a.data(h, t) = a.data(h, t) - a.data(col, t) * tmp
    '                Next
    '                '  x(h) = x(h) - x(col) * tmp
    '            End If
    '        Next
    '    Next
    '    For t = n - 1 To 0 Step -1
    '        If indfi(t) <> indco(t) Then
    '            For k = 0 To n - 1
    '                tmp = a.data(k, indfi(t)) = a.data(k, indco(t))
    '                a.data(k, indco(t)) = tmp
    '            Next
    '        End If
    '    Next
    'End Sub
    'Public Function multiply(ByVal mat1 As twoDMatrixLibrary.Matrix2Dsingle, ByVal mat2 As twoDMatrixLibrary.Matrix2Dsingle) As twoDMatrixLibrary.Matrix2Dsingle

    '    'we have 3 counters
    '    'first counters is the # of rows in the first twoDMatrixLibrary.Matrix
    '    'second counter is the # of cols in the first twoDMatrixLibrary.Matrix = #of rows in the 2nd ones
    '    'third counter is the # of cols in the 2nd twoDMatrixLibrary.Matrix

    '    If mat1.y <> mat2.x Then
    '        Return Nothing
    '    Else
    '        Dim ans As New twoDMatrixLibrary.Matrix2Dsingle(mat1.x, mat2.y)

    '        Dim a, b, c As Integer
    '        For a = 0 To mat1.x - 1
    '            For b = 0 To mat2.y - 1
    '                'a,b value is a row into b column
    '                ans.data(a, b) = 0
    '                For c = 0 To mat1.y - 1
    '                    ans.data(a, b) = ans.data(a, b) + (mat1.data(a, c) * mat2.data(c, b))
    '                Next
    '            Next
    '        Next
    '        Return ans
    '    End If
    'End Function
    'Public Function subtract(ByVal mat1 As twoDMatrixLibrary.Matrix2Dsingle, ByVal mat2 As twoDMatrixLibrary.Matrix2Dsingle) As twoDMatrixLibrary.Matrix2Dsingle
    '    Dim ans As New twoDMatrixLibrary.Matrix2Dsingle(mat1.x, mat1.y)
    '    Dim a, b As Integer
    '    For a = 0 To mat1.x - 1
    '        For b = 0 To mat2.y - 1
    '            ans.data(a, b) = mat1.data(a, b) - mat2.data(a, b)
    '        Next
    '    Next
    '    Return ans

    'End Function
    'Sub dialognalize() 'using pca

    '    'initialize eigenvector and eigenvalue matrices
    '    Me.EigenValue1 = New twoDMatrixLibrary.Matrix3DSingle(Nx, Ny, Nz)
    '    Me.EigenValue2 = New twoDMatrixLibrary.Matrix3DSingle(Nx, Ny, Nz)
    '    Me.EigenValue3 = New twoDMatrixLibrary.Matrix3DSingle(Nx, Ny, Nz)

    '    Me.EigenVector1 = New twoDMatrixLibrary.Matrix4Ddouble(Nx, Ny, Nz, 3)
    '    Me.EigenVector2 = New twoDMatrixLibrary.Matrix4Ddouble(Nx, Ny, Nz, 3)
    '    Me.EigenVector3 = New twoDMatrixLibrary.Matrix4Ddouble(Nx, Ny, Nz, 3)

    '    Dim a, b, c As Integer

    '    For a = 0 To Nx - 1
    '        For b = 0 To Ny - 1
    '            For c = 0 To Nz - 1

    '                Dim dtemp As New Mapack.Matrix(3, 3, 0)
    '                dtemp(0, 0) = Me.DiffusionMatrix.data(a, b, c, 0)
    '                dtemp(0, 1) = Me.DiffusionMatrix.data(a, b, c, 3)
    '                dtemp(0, 2) = Me.DiffusionMatrix.data(a, b, c, 5)
    '                dtemp(1, 0) = Me.DiffusionMatrix.data(a, b, c, 3)
    '                dtemp(1, 1) = Me.DiffusionMatrix.data(a, b, c, 1)
    '                dtemp(1, 2) = Me.DiffusionMatrix.data(a, b, c, 4)
    '                dtemp(2, 0) = Me.DiffusionMatrix.data(a, b, c, 5)
    '                dtemp(2, 1) = Me.DiffusionMatrix.data(a, b, c, 4)
    '                dtemp(2, 2) = Me.DiffusionMatrix.data(a, b, c, 2)

    '                Dim dtemp2 As New Mapack.EigenvalueDecomposition(dtemp)
    '                EigenValue1.data(a, b, c) = dtemp2.RealEigenvalues(0)
    '                EigenValue2.data(a, b, c) = dtemp2.RealEigenvalues(1)
    '                EigenValue1.data(a, b, c) = dtemp2.RealEigenvalues(2)
    '                EigenVector1.data(a, b, c, 0) = dtemp2.EigenvectorMatrix(0, 0)
    '                EigenVector1.data(a, b, c, 1) = dtemp2.EigenvectorMatrix(1, 0)
    '                EigenVector1.data(a, b, c, 2) = dtemp2.EigenvectorMatrix(2, 0)
    '                EigenVector2.data(a, b, c, 0) = dtemp2.EigenvectorMatrix(0, 1)
    '                EigenVector2.data(a, b, c, 1) = dtemp2.EigenvectorMatrix(1, 1)
    '                EigenVector2.data(a, b, c, 2) = dtemp2.EigenvectorMatrix(2, 1)
    '                EigenVector3.data(a, b, c, 0) = dtemp2.EigenvectorMatrix(0, 2)
    '                EigenVector3.data(a, b, c, 1) = dtemp2.EigenvectorMatrix(1, 2)
    '                EigenVector3.data(a, b, c, 2) = dtemp2.EigenvectorMatrix(2, 2)


    '            Next
    '        Next
    '    Next
    '    If MsgBox("Calculation of Eigenvectors and Eigenvalues complete. Press OK to start calculation of FA", MsgBoxStyle.OkCancel) = MsgBoxResult.Ok Then
    '        calculate_FA()
    '    End If

    'End Sub
    'Sub calculate_FA()
    '    Dim threshold As Single = 0.14
    '    Dim a, b, c As Integer
    '    Dim a1, a2, a3 As Single

    '    'initialize FA twoDMatrixLibrary.Matrix
    '    Me.FAMatrix = New twoDMatrixLibrary.Matrix3DSingle(Nx, Ny, Nz)

    '    For a = 0 To Nx - 1
    '        For b = 0 To Ny - 1
    '            For c = 0 To Nz - 1
    '                a1 = EigenValue1.data(a, b, c) 'largest
    '                a2 = EigenValue2.data(a, b, c)
    '                a3 = EigenValue3.data(a, b, c)
    '                Me.FAMatrix.data(a, b, c) = Math.Sqrt(((a1 - a2) ^ 2 + (a2 - a3) ^ 2 + (a3 - a1) ^ 2) / (2 * (a1 ^ 2 + a2 ^ 2 + a3 ^ 2)))

    '            Next
    '        Next
    '    Next


    '    MsgBox("Calculation of Fractional Anisotropy complete!")

    'End Sub
#End Region

#Region "Computation"
    Public Sub ComputeDiffusionMatrix()
        Dim x As New twoDMatrixLibrary.Matrix(Me.NumberOfGradients, 6)
        Dim i As Integer
        For i = 0 To Me.NumberOfGradients - 1
            x(i, 0) = Me.LoadedDTIDirections(i, 0) * Me.LoadedDTIDirections(i, 0)
            x(i, 1) = Me.LoadedDTIDirections(i, 1) * Me.LoadedDTIDirections(i, 1)
            x(i, 2) = Me.LoadedDTIDirections(i, 2) * Me.LoadedDTIDirections(i, 2)
            x(i, 3) = 2 * Me.LoadedDTIDirections(i, 0) * Me.LoadedDTIDirections(i, 1)
            x(i, 4) = 2 * Me.LoadedDTIDirections(i, 1) * Me.LoadedDTIDirections(i, 2)
            x(i, 5) = 2 * Me.LoadedDTIDirections(i, 0) * Me.LoadedDTIDirections(i, 2)
        Next
        Dim ls As twoDMatrixLibrary.Matrix = ((x.Transpose * x).Inverse) * x.Transpose
        Me.DiffusionMatrix = New Matrix4Ddouble(Me.Nx, Me.Ny, Me.Nz, 6)
        Dim res As New Matrix3DSingle(Nx, Ny, Nz)
        Dim trd As New Matrix3DSingle(Nx, Ny, Nz)
        Dim j, k As Integer
        For i = 0 To Nx - 1
            For j = 0 To Ny - 1
                For k = 0 To Nz - 1
                    Dim temp As New twoDMatrixLibrary.Matrix(Me.NumberOfGradients, 1)

                    Dim n As Integer
                    For n = 0 To Me.NumberOfGradients - 1
                        temp(n, 0) = Me.SubtractedMatrix.data(i, j, k, n)
                    Next

                    Dim dtemp As twoDMatrixLibrary.Matrix = (ls * temp)
                    Dim dif As twoDMatrixLibrary.Matrix = temp - (x * dtemp)

                    res.data(i, j, k) = CSng((Math.Sqrt((dif.Transpose * dif)(0, 0))) / ((Me.NumberOfGradients - 5) ^ 0.5))
                    trd.data(i, j, k) = CSng(dtemp(0, 0) + dtemp(1, 0) + dtemp(2, 0))

                    If trd.data(i, j, k) / res.data(i, j, k) > 0 Then 'check for threshold
                        DiffusionMatrix.data(i, j, k, 0) = dtemp(0, 0)
                        DiffusionMatrix.data(i, j, k, 1) = dtemp(1, 0)
                        DiffusionMatrix.data(i, j, k, 2) = dtemp(2, 0)
                        DiffusionMatrix.data(i, j, k, 3) = dtemp(3, 0)
                        DiffusionMatrix.data(i, j, k, 4) = dtemp(4, 0)
                        DiffusionMatrix.data(i, j, k, 5) = dtemp(5, 0)
                    End If
                Next
            Next
        Next
    End Sub

    Public Function ReturnComputedDiffusionmatrix() As Matrix4Ddouble
        Return Me.DiffusionMatrix
    End Function
    Public Sub DiagonalizeComputedDiffusionMatrix()
        Me.EigenValue1 = New Matrix3DSingle(Nx, Ny, Nz)
        Me.EigenValue2 = New Matrix3DSingle(Nx, Ny, Nz)
        Me.EigenValue3 = New Matrix3DSingle(Nx, Ny, Nz)
        Me.EigenVector1 = New Matrix4DdoubleS(Nx, Ny, Nz, 3)
        Me.EigenVector2 = New Matrix4DdoubleS(Nx, Ny, Nz, 3)
        Me.EigenVector3 = New Matrix4DdoubleS(Nx, Ny, Nz, 3)
        Dim Max As Double = 0
        Dim a, b, c As Integer

        For a = 0 To Nx - 1
            For b = 0 To Ny - 1
                For c = 0 To Nz - 1
                    Dim dtemp As New twoDMatrixLibrary.Matrix(3, 3, 0)
                    dtemp(0, 0) = DiffusionMatrix.data(a, b, c, 0)
                    dtemp(0, 1) = DiffusionMatrix.data(a, b, c, 3)
                    dtemp(0, 2) = DiffusionMatrix.data(a, b, c, 5)
                    dtemp(1, 0) = DiffusionMatrix.data(a, b, c, 3)
                    dtemp(1, 1) = DiffusionMatrix.data(a, b, c, 1)
                    dtemp(1, 2) = DiffusionMatrix.data(a, b, c, 4)
                    dtemp(2, 0) = DiffusionMatrix.data(a, b, c, 5)
                    dtemp(2, 1) = DiffusionMatrix.data(a, b, c, 4)
                    dtemp(2, 2) = DiffusionMatrix.data(a, b, c, 2)
                    Dim dtemp2 As New EigenvalueDecomposition(dtemp)
                    EigenValue1.data(a, b, c) = CSng(dtemp2.RealEigenvalues(0))
                    EigenValue2.data(a, b, c) = CSng(dtemp2.RealEigenvalues(1))
                    EigenValue3.data(a, b, c) = CSng(dtemp2.RealEigenvalues(2))
                    If EigenValue1(a, b, c) > Max Then Max = EigenValue1(a, b, c)
                    EigenVector1.data(a, b, c, 0) = dtemp2.EigenvectorMatrix(0, 0)
                    EigenVector1.data(a, b, c, 1) = dtemp2.EigenvectorMatrix(1, 0)
                    EigenVector1.data(a, b, c, 2) = dtemp2.EigenvectorMatrix(2, 0)
                    EigenVector2.data(a, b, c, 0) = dtemp2.EigenvectorMatrix(0, 1)
                    EigenVector2.data(a, b, c, 1) = dtemp2.EigenvectorMatrix(1, 1)
                    EigenVector2.data(a, b, c, 2) = dtemp2.EigenvectorMatrix(2, 1)
                    EigenVector3.data(a, b, c, 0) = dtemp2.EigenvectorMatrix(0, 2)
                    EigenVector3.data(a, b, c, 1) = dtemp2.EigenvectorMatrix(1, 2)
                    EigenVector3.data(a, b, c, 2) = dtemp2.EigenvectorMatrix(2, 2)

                Next
            Next
        Next


    End Sub
    Function ReturnEigenVector(ByVal EigenVectorCount_1_2_or_3 As Integer) As Matrix4DdoubleS
        Select Case EigenVectorCount_1_2_or_3
            Case 1
                Return Me.EigenVector1
            Case 2
                Return Me.EigenVector2
            Case 3
                Return Me.EigenVector3
            Case Else
                Return Nothing
        End Select
    End Function
    Function ReturnEigenVector(ByVal EigenVectorCount_1_2_or_3 As Integer, ByVal x As Integer, ByVal y As Integer, ByVal z As Integer, ByVal count As Integer) As Integer
        Select Case EigenVectorCount_1_2_or_3
            Case 1
                Return CInt(Me.EigenVector1.data(x, y, z, count))
            Case 2
                Return CInt(Me.EigenVector2.data(x, y, z, count))
            Case 3
                Return CInt(Me.EigenVector3.data(x, y, z, count))
            Case Else
                Return Nothing
        End Select
    End Function
    Function ReturnEigenValue(ByVal EigenValueCount_1_2_or_3 As Integer) As Matrix3DSingle
        Select Case EigenValueCount_1_2_or_3
            Case 1
                Return Me.EigenValue1
            Case 2
                Return Me.EigenValue2
            Case 3
                Return Me.EigenValue3
            Case Else
                MessageBox.Show("Parameter: " & EigenValueCount_1_2_or_3)
                Return Nothing
        End Select
    End Function
    Function ReturnEigenValue(ByVal EigenValueCount_1_2_or_3 As Integer, ByVal x As Integer, ByVal y As Integer, ByVal z As Integer) As Integer
        Select Case EigenValueCount_1_2_or_3
            Case 1
                Return CInt(Me.EigenValue1.data(x, y, z))
            Case 2
                Return CInt(Me.EigenValue2.data(x, y, z))
            Case 3
                Return CInt(Me.EigenValue3.data(x, y, z))
            Case Else
                Return Nothing
        End Select
    End Function
    Sub CalculateFA()
        Dim a, b, c As Integer
        Dim a1, a2, a3, a4 As Double

        'initialize FA twoDMatrixLibrary.Matrix
        Dim FAmax As Single = 0
        Me.FAMatrix = New Matrix3DSingle(Nx, Ny, Nz)
        For a = 0 To Nx - 1
            For b = 0 To Ny - 1
                For c = 0 To Nz - 1
                    a1 = EigenValue1.data(a, b, c)
                    a2 = EigenValue2.data(a, b, c)
                    a3 = EigenValue3.data(a, b, c)
                    a4 = (a1 + a2 + a3) / 3
                    If (a4 <> 0) Then
                        Me.FAMatrix(a, b, c) = CSng(3 * (((a1 - a4) ^ 2) + ((a2 - a4) ^ 2) + ((a3 - a4) ^ 2)))
                        Me.FAMatrix(a, b, c) = CSng(Me.FAMatrix(a, b, c) / (2 * ((a1 ^ 2) + (a2 ^ 2) + (a3 ^ 2))))
                        Me.FAMatrix(a, b, c) = CSng(Me.FAMatrix(a, b, c) ^ 0.5)
                    Else
                        Me.FAMatrix(a, b, c) = 0
                    End If

                    If Me.FAMatrix(a, b, c) > 1 Then Me.FAMatrix(a, b, c) = 1
                Next
            Next
        Next
        'For a = 0 To Nx - 1
        '    For b = 0 To Ny - 1
        '        For c = 0 To Nz - 1
        '            FAMatrix(a, b, c) /= FAmax
        '        Next
        '    Next
        'Next
        RaiseEvent FACalculationComplete(Nothing, EventArgs.Empty)
        '  Dim ass As New ParRecLoader.DataScaler
        '  ass.Scale3DData(FAMatrix, 0, 256)
    End Sub

    Function ReturnFAMatrixSingle() As Matrix3DSingle
        Return DTIComputer.FA
    End Function
    Function ReturnFAMatrixInt16() As Matrix3DInt16
        Dim a, b, c As Integer
        Dim ans As New Matrix3DInt16(Nx, Ny, Nz)
        For a = 0 To Nx - 1
            For b = 0 To Ny - 1
                For c = 0 To Nz - 1
                    ans.data(a, b, c) = Convert.ToInt16(Me.FAMatrix.data(a, b, c))
                Next
            Next
        Next
        Return ans
    End Function
    Function ReturnFAValue(ByVal x As Integer, ByVal y As Integer, ByVal z As Integer) As Single
        Return Me.FAMatrix.data(x, y, z)
    End Function
#End Region

#Region "constructors"
    Sub New()
        InitializeComponent()
    End Sub

    Sub New(ByVal RawData As Matrix4DInt16, ByVal x As Integer, ByVal y As Integer, ByVal z As Integer)
        ' This call is required by the Windows Form Designer.
        InitializeComponent()
        ' Add any initialization after the InitializeComponent() call.
        Me.RawDTIData = New Matrix4DInt16(RawData.data)
        Nx = Me.RawDTIData.i
        Ny = Me.RawDTIData.j
        Nz = Me.RawDTIData.k
    End Sub
#End Region

#Region "Open Save Directions DIalog"
    Public Function SeperateNumbersFromString(ByVal a_s As String) As String()
        Dim a() As Char
        a = a_s.ToCharArray
        Dim i, j As Integer
        i = 0
        Dim dir(2) As String
        While a(i) <> " "
            i = i + 1
        End While
        dir(0) = Mid(a_s, 1, i) 'get the first direction
        i = i + 1               'for " "
        j = i
        While a(j) <> " "
            j = j + 1
        End While
        dir(1) = Mid(a_s, i, j - i + 1)
        dir(2) = Mid(a_s, j + 1, Len(a_s) - j)
        Return dir
    End Function
    Private Sub ClearButtonClick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btn_clear.Click
        txt_gradient_strength.Text = "0"
        txt_small_delta.Text = "0".ToString
        txt_big_delta.Text = "0"
        dir_list.Items.Clear()
        dir_counter.Value = 0
    End Sub
    Private Sub AddDirectionEntry(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles txt_add_dir.Click
        Dim a_s As String
        a_s = InputBox.Show("Enter new direction (xx.xx xx.xx xx.xx):", "Add direction").Text
        Dim dir(2) As String
        dir = SeperateNumbersFromString(a_s)

        If dir(0) = "" Or dir(0) = " " Or dir(1) = "" Or dir(1) = " " Or dir(2) = "" Or dir(2) = " " Then
            MsgBox("Improper format entered!", MsgBoxStyle.OkOnly)
        Else
            dir_counter.UpButton()
            Dim item As ListViewItem
            item = New ListViewItem(dir)
            dir_list.Items.Add(item)
        End If
    End Sub
    Private Sub DeleteDirectionEntry(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles txt_del_dir.Click
        If dir_list.SelectedItems.Count <> 1 Then
            MsgBox("Please select 1 direction to delete!", MsgBoxStyle.OkOnly)
        Else
            dir_counter.DownButton()
            dir_list.Items(dir_list.SelectedItems(0).Index).Remove()
        End If
    End Sub
    Private Sub LoadFromXMLButtonCLick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btn_load_from_xml.Click
        OpenFileDialog1.Filter = "XML Brainfibercraft File (*.xml)|*.xml|All files (*.*)|*.*"
        OpenFileDialog1.Title = "Choose the XML File"
        OpenFileDialog1.ShowDialog()
    End Sub
    Private Sub OpenFileDialog1_FileOk(ByVal sender As System.Object, ByVal e As System.ComponentModel.CancelEventArgs) Handles OpenFileDialog1.FileOk
        Try
            ScannerParamFile(OpenFileDialog1.FileName)
        Catch ex As Exception
            MsgBox("File does not exist or Invalid File!")
        End Try
    End Sub

    Public Function ScannerParamFile(ByVal filename As String) As XmlDocument
        Dim xmlDoc As XmlDocument
        xmlDoc = New XmlDocument
        xmlDoc.Load(filename)

        Dim grad_str, big_d, small_d, dir_count As XmlNode
        If xmlDoc.DocumentElement.HasChildNodes Then

            'first 
            grad_str = xmlDoc.DocumentElement.FirstChild
            If grad_str.Name = "gradient_strength" Then
                txt_gradient_strength.Text = grad_str.ChildNodes(0).Value
            End If

            big_d = xmlDoc.ChildNodes(1).ChildNodes(1)
            If big_d.Name = "big_d" Then
                txt_big_delta.Text = big_d.ChildNodes(0).Value
            End If
            'next
            small_d = xmlDoc.ChildNodes(1).ChildNodes(2)
            If small_d.Name = "small_d" Then
                txt_small_delta.Text = small_d.ChildNodes(0).Value
            End If

            dir_count = xmlDoc.ChildNodes(1).ChildNodes(3)
            Dim dir_c As Integer = CInt(CDbl(dir_count.ChildNodes(0).Value) * 3)

            Dim i As Integer
            For i = 1 To dir_c Step 3
                Dim dir(2) As String
                Dim dir_node As XmlNode
                dir_node = xmlDoc.ChildNodes(1).ChildNodes(3 + i)
                Dim direction As String = dir_node.ChildNodes(0).Value
                dir(0) = direction
                dir_node = xmlDoc.ChildNodes(1).ChildNodes(i + 4)
                direction = dir_node.ChildNodes(0).Value
                dir(1) = direction
                dir_node = xmlDoc.ChildNodes(1).ChildNodes(i + 5)
                direction = dir_node.ChildNodes(0).Value
                dir(2) = direction
                Dim item As ListViewItem
                item = New ListViewItem(dir)
                dir_list.Items.Add(item)
                dir_counter.UpButton()
            Next i
        End If
        Return xmlDoc
    End Function
    Private Sub SaveToXMLButtonClick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btn_save_to_xml.Click
        Dim xmlDoc As XmlDocument
        xmlDoc = New XmlDocument

        'main document node to be appended into
        Dim parentnode As XmlElement


        'new nodes
        Dim childnode(3) As XmlElement


        parentnode = xmlDoc.DocumentElement
        parentnode = xmlDoc.CreateElement("parameters")

        Dim start As XmlDeclaration
        start = xmlDoc.CreateXmlDeclaration("1.0", "UTF-8", "yes")
        xmlDoc.AppendChild(start)

        childnode(0) = xmlDoc.CreateElement("gradient_strength")
        childnode(1) = xmlDoc.CreateElement("big_d")
        childnode(2) = xmlDoc.CreateElement("small_d")
        childnode(3) = xmlDoc.CreateElement("dir_count")

        Dim xmlgd, xmlbd, xmlsd, xmldc As XmlText

        xmlgd = xmlDoc.CreateTextNode("gradient_strength")
        xmlgd.Value = txt_gradient_strength.Text
        childnode(0).AppendChild(xmlgd)
        parentnode.AppendChild(childnode(0))

        xmlbd = xmlDoc.CreateTextNode("big_d")
        xmlbd.Value = txt_big_delta.Text
        childnode(1).AppendChild(xmlbd)
        parentnode.AppendChild(childnode(1))

        xmlsd = xmlDoc.CreateTextNode("small_d")
        xmlsd.Value = txt_small_delta.Text
        childnode(2).AppendChild(xmlsd)
        parentnode.AppendChild(childnode(2))

        xmldc = xmlDoc.CreateTextNode("dir_count")
        xmldc.Value = dir_counter.Value.ToString
        childnode(3).AppendChild(xmldc)
        parentnode.AppendChild(childnode(3))

        Dim i As Integer
        If dir_counter.Value > 0 Then
            For i = 1 To CInt(dir_counter.Value)
                Dim xmlvd1, xmlvd2, xmlvd3 As XmlText
                Dim childnodev(2) As XmlElement
                Dim item As ListViewItem
                item = dir_list.Items(i - 1)
                childnodev(0) = xmlDoc.CreateElement("value1")
                childnodev(1) = xmlDoc.CreateElement("value2")
                childnodev(2) = xmlDoc.CreateElement("value3")
                xmlvd1 = xmlDoc.CreateTextNode("value1")
                xmlvd1.Value = item.SubItems(0).Text
                childnodev(0).AppendChild(xmlvd1)
                parentnode.AppendChild(childnodev(0))
                xmlvd2 = xmlDoc.CreateTextNode("value2")
                xmlvd2.Value = item.SubItems(0).Text
                childnodev(1).AppendChild(xmlvd2)
                parentnode.AppendChild(childnodev(1))
                xmlvd3 = xmlDoc.CreateTextNode("value3")
                xmlvd3.Value = item.SubItems(0).Text
                childnodev(2).AppendChild(xmlvd3)
                parentnode.AppendChild(childnodev(2))
            Next

            'add parentnode to main document parent node and save
            xmlDoc.AppendChild(parentnode)
            xmlDoc.Save("c:\DTI Parameters.xml")
        Else
            MsgBox("Direction List is empty")
        End If
    End Sub
#End Region

    Private Sub Button2_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button2.Click
        Me.Close()
    End Sub

    Dim DTIComputer As ComputeDTI

    Public Function ReturnFA() As Matrix3DSingle
        If choice Then
            Return DTIComputer.FA
        Else
            Return Me.FAMatrix
        End If
    End Function
    Public Function ReturnAD() As Matrix3DSingle
        Return DTIComputer.AD
    End Function
    Public Function ReturnLinearlization() As Matrix3DSingle
        Return DTIComputer.Linearization
    End Function
    Public Function ReturnEigenValue(ByVal Number As ComputeDTI.Order) As Matrix3DSingle
        If choice Then
            Select Case Number
                Case ComputeDTI.Order.Big
                    Return DTIComputer.LambdaBig
                Case ComputeDTI.Order.Medium
                    Return DTIComputer.LambdaMed
                Case ComputeDTI.Order.Small
                    Return DTIComputer.LambdaSmall
            End Select
            Return Nothing
        Else
            Select Case Number
                Case ComputeDTI.Order.Big
                    Return Me.EigenValue1
                Case ComputeDTI.Order.Medium
                    Return Me.EigenValue2
                Case ComputeDTI.Order.Small
                    Return Me.EigenValue3
            End Select
            Return Nothing
        End If
    End Function
    Public Function ReturnEigenVector(ByVal Number As ComputeDTI.Order) As Matrix4DdoubleS
        If choice Then
            Select Case Number
                Case ComputeDTI.Order.Big
                    Return DTIComputer.VectorBig
                Case ComputeDTI.Order.Medium
                    Return DTIComputer.VectorMed
                Case ComputeDTI.Order.Small
                    Return DTIComputer.VectorSmall
            End Select
            Return Nothing
        Else
            Select Case Number
                Case ComputeDTI.Order.Big
                    Return Me.EigenVector1
                Case ComputeDTI.Order.Medium
                    Return Me.EigenVector2
                Case ComputeDTI.Order.Small
                    Return Me.EigenVector3
            End Select
            Return Nothing
        End If
    End Function
    Private choice As Boolean = False
    Private Sub OKButtonClick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        Me.NumberOfGradients = CInt(dir_counter.Value)
        Me.GradientStrength = CSng(Mid(txt_gradient_strength.Text.Trim, 1))
        Me.SmallDelta = CSng(Mid(txt_small_delta.Text.Trim, 1))
        Me.BigDelta = CSng(Mid(txt_big_delta.Text.Trim, 1))
        If Me.GradientStrength <> 0.0F And Me.SmallDelta <> 0.0F And BigDelta <> 0.0F And Me.NumberOfGradients > 6 Then
            LoadedDTIDirections = New twoDMatrixLibrary.Matrix(NumberOfGradients, 3)
            Dim i As Integer
            For i = 0 To NumberOfGradients - 1
                Me.LoadedDTIDirections(i, 0) = CSng(dir_list.Items(i).SubItems(0).Text.Trim)
                Me.LoadedDTIDirections(i, 1) = CSng(dir_list.Items(i).SubItems(1).Text.Trim)
                Me.LoadedDTIDirections(i, 2) = CSng(dir_list.Items(i).SubItems(2).Text.Trim)
            Next
            Me.Opacity = 0.5
        End If
        Dim temp As Boolean = False
        If MsgBox("Directions Loaded! Closing window and proceeding with FA calculation." & vbNewLine & _
            "No of Gradients entered: " & Me.NumberOfGradients & ". Number of Gradients Detected: " & Me.RawDTIData.l & vbNewLine & _
                  "Do you want to save temporary matrices/files for debugging (Yes/No) ?", MsgBoxStyle.YesNo, "Save Debugging Information") = MsgBoxResult.Yes Then temp = True

        Me.Visible = False
        If MsgBox("Use 6-directions only to calculate FA?", MsgBoxStyle.YesNo) = MsgBoxResult.No Then
            choice = True
            DTIComputer = New ComputeDTI(Me.BigDelta, Me.SmallDelta, Me.GradientStrength, Me.NumberOfGradients, Me.LoadedDTIDirections, Me.RawDTIData, temp)
            DTIComputer.Go()
            RaiseEvent FACalculationComplete(Nothing, EventArgs.Empty)
        Else
            choice = False
            Me.StartDTIAnalysis()
            Me.ComputeDiffusionMatrix()
            Me.DiagonalizeComputedDiffusionMatrix()
            Me.CalculateFA()
            RaiseEvent FACalculationComplete(Nothing, EventArgs.Empty)
        End If
    End Sub

End Class


#Region "New DTI Computation Class based on Bear"

Public Class ComputeDTI
    '1. Compute B Matrix
    '2. Create Brain Mask from A0 data and apply it to everything
    '3. Take Ratio of each direction with A0 data
    '4. Take log of ratio in each slice where ratio is not 0.
    '5. Diagnoalize the diffusion tensor for each voxel
    '6. Reconstruct V1,V2,V3,AD,FA,L1,L2,L3,DD

#Region "Constructors & Initialization"
    Dim SmallDelta, BigDelta, GradientStrength As Double
    Dim NumberofGradients As Integer
    Public GradientDirections, Bmatrix As twoDMatrixLibrary.Matrix
    Public RawData As Matrix4DInt16
    Dim A0Data As Matrix3DInt16
    Public Enum Order
        Big = 2
        Medium = 1
        Small = 0
    End Enum
    Public LambdaBig, LambdaMed, LambdaSmall, FA, AD, Linearization As Matrix3DSingle
    Public VectorBig, VectorMed, VectorSmall As Matrix4DdoubleS

    Dim SaveAllTempMatrices As Boolean = False

    Public Sub New(ByVal BigDelta As Double, ByVal SmallDelta As Double, ByVal GradientStrength As Double, ByVal NumberofGradients As Integer, ByVal GradientDirections As Matrix, ByVal Data As Matrix4DInt16, ByVal SaveAllTempMatrices As Boolean)
        Me.BigDelta = BigDelta
        Me.SmallDelta = SmallDelta
        Me.GradientStrength = GradientStrength
        Me.GradientDirections = GradientDirections
        Me.NumberofGradients = NumberofGradients
        Me.SaveAllTempMatrices = SaveAllTempMatrices
        Me.RawData = Data
    End Sub

#End Region

#Region "1. Compute B Matrix"
    Sub ComputeBMatrix()
        Dim AlphaSquare As Double = (26750.0 * 26750.0) * (GradientStrength * GradientStrength) * (BigDelta * BigDelta) * (BigDelta - SmallDelta / 3)
        AlphaSquare = Double.Parse(InputBox.Show("Enter a different b-value", "B-value", AlphaSquare.ToString()).Text)
        Dim Temp As Matrix3DSingle = New Matrix3DSingle(3, 3, Me.NumberofGradients)
        For n As Integer = 0 To Me.NumberofGradients - 1
            Dim nslice As twoDMatrixLibrary.Matrix = Me.GradientDirections.Submatrix(n, n, 0, 2).Transpose * Me.GradientDirections.Submatrix(n, n, 0, 2)
            For i As Integer = 0 To 2
                For j As Integer = 0 To 2
                    Temp(i, j, n) = CSng(nslice(i, j))
                Next
            Next
        Next
        Me.Bmatrix = New Matrix(Me.NumberofGradients, 6)
        For i As Integer = 0 To Me.NumberofGradients - 1
            For j As Integer = 0 To 2
                Bmatrix(i, j) = -1 * Temp(j, j, i) * AlphaSquare
            Next
            Bmatrix(i, 3) = -2 * Temp(0, 1, i) * AlphaSquare
            Bmatrix(i, 4) = -2 * Temp(0, 2, i) * AlphaSquare
            Bmatrix(i, 5) = -2 * Temp(1, 2, i) * AlphaSquare
        Next
        If Me.SaveAllTempMatrices Then Bmatrix.SaveAs(System.IO.Path.GetTempPath & "\bmatrix.txt")
        Temp.Dispose()
    End Sub
#End Region
    'rest 3 steps are performed slicewise in 2D
#Region "4. Take log of ratio in each slice where ratio is not 0. 5. Diagnoalize the diffusion tensor for each voxel 6. Reconstruct V1,V2,V3,AD,FA,L1,L2,L3,DD"
    Sub PerformDTComputation()
        GC.Collect(1)
        GC.WaitForPendingFinalizers()

        Dim DiffusionWeightedImageset As Matrix3DSingle = New Matrix3DSingle(Me.RawData.i, Me.RawData.j, Me.NumberofGradients)
        Dim AA As Matrix = New Matrix(Me.NumberofGradients, 1)

        Dim BI As Matrix = (Me.Bmatrix.Transpose * Me.Bmatrix).Inverse * Bmatrix.Transpose
        If Me.SaveAllTempMatrices Then BI.SaveAs(System.IO.Path.GetTempPath & "\BI.txt")

        Dim temp As Matrix, Stdev As Single
        Dim DiffusionTensorOfVoxel As Matrix = New Matrix(3, 3)
        Dim Diagonalizer As twoDMatrixLibrary.EigenvalueDecomposition
        Me.LambdaBig = New Matrix3DSingle(Me.RawData.i, Me.RawData.j, Me.RawData.k)
        Me.LambdaMed = New Matrix3DSingle(Me.RawData.i, Me.RawData.j, Me.RawData.k)
        Me.LambdaSmall = New Matrix3DSingle(Me.RawData.i, Me.RawData.j, Me.RawData.k)
        Me.FA = New Matrix3DSingle(Me.RawData.i, Me.RawData.j, Me.RawData.k)
        Me.AD = New Matrix3DSingle(Me.RawData.i, Me.RawData.j, Me.RawData.k)
        Me.Linearization = New Matrix3DSingle(Me.RawData.i, Me.RawData.j, Me.RawData.k)
        Me.VectorBig = New Matrix4DdoubleS(Me.RawData.i, Me.RawData.j, Me.RawData.k, 3)
        Me.VectorMed = New Matrix4DdoubleS(Me.RawData.i, Me.RawData.j, Me.RawData.k, 3)
        Me.VectorSmall = New Matrix4DdoubleS(Me.RawData.i, Me.RawData.j, Me.RawData.k, 3)
        Dim counter As Integer = 0
        Dim map1 As Boolean = (MsgBox("Do you want to map all values greater than 1 to 1?", MsgBoxStyle.YesNo, "") = MsgBoxResult.Yes)

        Dim b_tag As Integer = 0
        b_tag = Integer.Parse(InputBox.Show("Which is the B=0 Matrix? [starts from 0]", "B=0", b_tag.ToString).Text)

        For k As Integer = 0 To Me.RawData.k - 1
            For i As Integer = 0 To Me.RawData.i - 1
                For j As Integer = 0 To Me.RawData.j - 1
                    ' for each slice
                    'Step 4
                    counter = 0
                    For l As Integer = 0 To Me.NumberofGradients - 1 'MAJOR change. 
                        Me.RawData.data(i, j, k, l) = CShort(IIf(Me.RawData.data(i, j, k, l) = 0, 1, Me.RawData.data(i, j, k, l)))
                        DiffusionWeightedImageset(i, j, l) = CSng((Me.RawData.data(i, j, k, l + 1)) / Me.RawData.data(i, j, k, 0))

                        If DiffusionWeightedImageset(i, j, l) <> 0 Then
                            AA(l, 0) = Math.Log(DiffusionWeightedImageset(i, j, l), Math.E)
                        Else
                            AA(l, 0) = 0
                            counter += 1
                        End If
                    Next

                    If (Me.NumberofGradients - counter) > 6 Then
                        'all values in AA are defined now
                        temp = BI * AA

                        DiffusionTensorOfVoxel(0, 0) = temp(0, 0)
                        DiffusionTensorOfVoxel(0, 1) = temp(3, 0)
                        DiffusionTensorOfVoxel(0, 2) = temp(4, 0)
                        DiffusionTensorOfVoxel(1, 0) = temp(3, 0)
                        DiffusionTensorOfVoxel(1, 1) = temp(1, 0)
                        DiffusionTensorOfVoxel(1, 2) = temp(5, 0)
                        DiffusionTensorOfVoxel(2, 0) = temp(4, 0)
                        DiffusionTensorOfVoxel(2, 1) = temp(5, 0)
                        DiffusionTensorOfVoxel(2, 2) = temp(2, 0)
                        If Me.SaveAllTempMatrices Then DiffusionTensorOfVoxel.SaveAs(System.IO.Path.GetTempPath & "voxel_" & i & "_" & j & "_" & k & ".txt")

                        'Step 5
                        Diagonalizer = New twoDMatrixLibrary.EigenvalueDecomposition(DiffusionTensorOfVoxel)

                        'Step 6
                        Me.LambdaBig(i, j, k) = CSng(Diagonalizer.RealEigenvalues(0))
                        Me.LambdaMed(i, j, k) = CSng(Diagonalizer.RealEigenvalues(1))
                        Me.LambdaSmall(i, j, k) = CSng(Diagonalizer.RealEigenvalues(2))

                        Me.AD(i, j, k) = (Me.LambdaBig(i, j, k) + Me.LambdaMed(i, j, k) + Me.LambdaSmall(i, j, k)) / 3

                        Stdev = CSng(((((Me.LambdaBig(i, j, k) - Me.AD(i, j, k)) ^ 2) + _
                                       ((Me.LambdaMed(i, j, k) - Me.AD(i, j, k)) ^ 2) + _
                                       ((Me.LambdaSmall(i, j, k) - Me.AD(i, j, k)) ^ 2)) / 3) ^ (0.5))


                        If LambdaBig(i, j, k) ^ 2 + LambdaMed(i, j, k) ^ 2 + LambdaSmall(i, j, k) ^ 2 = 0 Then
                            FA(i, j, k) = 0
                        Else
                            FA(i, j, k) = CSng(3 * Stdev / ((2 * ((LambdaBig(i, j, k) ^ 2) + (LambdaMed(i, j, k) ^ 2) + (LambdaSmall(i, j, k) ^ 2))) ^ 0.5))
                            If map1 And FA(i, j, k) > 1 Then FA(i, j, k) = 1
                            Me.Linearization(i, j, k) = (LambdaBig(i, j, k) - LambdaMed(i, j, k)) * 100 / (3 * AD(i, j, k))

                        End If

                        'eigenvectors now
                        Me.VectorBig.data(i, j, k, 0) = Diagonalizer.EigenvectorMatrix(0, 0)
                        Me.VectorBig.data(i, j, k, 1) = Diagonalizer.EigenvectorMatrix(1, 0)
                        Me.VectorBig.data(i, j, k, 2) = Diagonalizer.EigenvectorMatrix(2, 0)
                        Me.VectorMed.data(i, j, k, 0) = Diagonalizer.EigenvectorMatrix(0, 1)
                        Me.VectorMed.data(i, j, k, 1) = Diagonalizer.EigenvectorMatrix(1, 1)
                        Me.VectorMed.data(i, j, k, 2) = Diagonalizer.EigenvectorMatrix(2, 1)
                        Me.VectorSmall.data(i, j, k, 0) = Diagonalizer.EigenvectorMatrix(0, 2)
                        Me.VectorSmall.data(i, j, k, 1) = Diagonalizer.EigenvectorMatrix(1, 2)
                        Me.VectorSmall.data(i, j, k, 2) = Diagonalizer.EigenvectorMatrix(2, 2)
                    End If
                Next
            Next
        Next
    End Sub
#End Region

    Public Sub Go()
        Me.ComputeBMatrix()
        Me.PerformDTComputation()
        If Me.SaveAllTempMatrices Then MsgBox("Temp files saved to: " & System.IO.Path.GetTempPath)
    End Sub
End Class

#End Region

